runPBS <- function(.output,
                   .priors,
                   .bsiter,
                   .rc = NULL,
                   .thresh = 1e-6,
                   .asEM = FALSE
                   ) {

    ## #############
    ## BS Containers
    ## #############

    aX <- array(data = NA, dim = c(.output$n, .output$d, .bsiter))
    aBeta <- array(data = NA, dim = c(.output$j, .output$d + 1, .bsiter))
    vDebug <- numeric(.bsiter)

    ## ####################
    ## Fitted Probabilities
    ## ####################

    x2 <- cbind(1,
                .output$means$x
                )
    mYSH <- x2 %*% t(.output$means$beta)
    mProbs <- pnorm(mYSH)


    for (bs in 1:.bsiter) {
        mYbs <- matrix(rbinom(prod(dim(mProbs)), 1, mProbs),
                       nrow = nrow(mProbs)
                       )

        rcbs <- rollcall(mYbs)
        newrcbs <- convertRC(rcbs)
        if (!is.null(.rc)) {
            isunobs <- newrcbs$votes %in% c(.rc$codes$missing, .rc$codes$notInLegis)
            newrcbs$votes[isunobs] <- 0
            ## table(newrcbs$votes)
            ## table(newrc$votes)
        }

        loutbs <- fastest(.rc = newrcbs,
                          .starts = (list(x = .output$means$x,
                                          beta = .output$means$beta[, 2, drop = FALSE],
                                          alpha = .output$means$beta[, 1, drop = FALSE]
                                          )
                                     ),
                          .priors = .priors,
                          .control = {list(
                              threads = 1,
                              verbose = TRUE,
                              thresh = .thresh,
                              maxit = 4000,
                              convtype = 1,
                              asEM = .asEM
                              )
                                  }
                          )

        vDebug[bs] <- cor(loutbs$means$x, .output$means$x)
        aX[, , bs] <- loutbs$means$x
        aBeta[, , bs] <- loutbs$means$beta
    }

    mXvar <- apply(aX, c(1, 2), var)
    mBetavar <- apply(aBeta, c(1, 2), var)

    ## #######
    ## Returns
    ## #######

    ret <- list(x = mXvar,
                beta = mBetavar,
                debug = vDebug
                )

    return(ret)
}
